rm(list = ls())
library(reshape2)
library(ggplot2)
library(magrittr)
library(stringr)
devtools::load_all()

# 参数设置
beta <- .98
delta <-  .1
theta <- .36


# 状态变量格点化
k0 <- k <- seq(4,6,0.06)
vlast <- rep(0.2,length(k))
picdata <- data.frame(k0 = k0, v = vlast)

#
maxfun <- function(cc,k,vlast){
  knext <- k^theta + (1-delta)*k - cc
  g <- signal::interp1(k0,vlast,knext,'linear')

  val <- log(cc)+beta*g

  return(-val)
}

c_opt <- v <- numeric(length(k))
for (i in 1:10) {
  for (j in 1:length(k)) {
    ans <- optimize(maxfun,k=k[j],interval = c(0.1,6.2),vlast = vlast)

    v[j] <- -ans$objective
  }

  # 每48次迭代存储一下
  if (i %% 1 == 0){
    print(i)
    picdata[,paste('v',as.character(i), sep = '')] <- v
    # picdata[,paste('k_opt',as.character(i), sep = '')] <- k_opt
  }
  vlast <- v
}



